43 
Oh 



,c 



Cosmological Perturbations at Second Order 

and 

Recombination Perturbed 

Leonardo Senatore a ' 6 ' c , Svetlin Tassev c and Matias Zaldarriaga 6 

a School of Natural Sciences, Institute for Advanced Study, 

00 Olden Lane, Princeton, NJ 08540, USA 

O 

b Jefferson Physical Laboratory, Harvard University, Cambridge, MA 02138, USA 

O 

Q 

00 



Center for Astrophysics, Harvard University, Cambridge, MA 02138, USA 



Abstract 

We derive the full set of second-order equations governing the evolution of cosmological perturba- 
tions, including the effects of the first-order electron number density perturbations, 5 e . We provide 
a detailed analysis of the perturbations to the recombination history of the universe and show that 
a perturbed version of the Peebles effective 3-level atom is sufficient for obtaining the evolution of S e 
i— i for comoving wavenumbers smaller than 1 Mpc . We calculate rigorously the perturbations to the 
Lya escape probability and show that to a good approximation it is governed by the local baryon 
velocity divergence. For modes shorter than the photon diffusion scale, we find that 5 e is enhanced 
during recombination by a factor of roughly 5 relative to other first-order quantities sourcing the 
^sO CMB anisotropies at second order. Using these results, in a companion paper we calculate the CMB 
bispectrum generated during recombination. 
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1 Introduction 

> 

^ In the last few years there has been great progress in understanding the non-Gaussianity of the 
primordial spectrum of density fluctuations. Starting from the first full computation of the non- 
Gaussian features in single field slow roll inflation [1] [2] , several alternative models have been pro- 
posed that produce a large and in principle detectable level of non-Gaussianities through different 
mechanisms for generating density fluctuations in the quasi de Sitter inflationary phase, both in the 
case of single field inflation j3j HJ EJ EJ [7J , and in the case of multi-field inflation [U [9]. Further, it 
has been at last developed a consistent model of bouncing universe |10[ [TTJ [12] which, though clearly 
less compelling than inflation, predicts a possibly detectable amount of non-Gaussianity [12]. If 
large non-Gaussianities will be detected, we would have to abandon the standard slow-roll inflation 
picture, but would be left with new information to understand the dynamics of the inflaton |13[ [7j. 

At the same time, from the experimental side, the WMAP satellite has allowed for a huge 
improvement in our measurement of the properties of the CMB. Observations seem to confirm the 
generic predictions of standard slow roll inflation |14j . Limits on the primordial non-Gaussianity 
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of the CMB have been significantly improved |15[ [T4] . but for the moment the data are consistent 
with a non-Gaussian signal. Recently it has been realized that even large scale structure surveys are 
relevantly sensitive to at least some particular kind of non-Gaussianities |16| [T7] , giving constraints 
comparable to those from WMAP |17j . 

Oversimplifying, the limits on non-Gaussianities are set on the parameter /nl such that the 
curvature £ that we observe in the CMB is given by: 



where C g {x) is a gaussian random variable Q So far we have a l-a error on /nl of order A/nl — 30 
from WMAP [14 a , and soon the Planck satellite is expected to improve this limit to A/nl ^ 4 [TH] . 
However, it is clear by looking at eq. ([TJ that the corrections due to the non-linearities of general 
relativity and of the Boltzmann equations that regulate the fluid dynamics are very naively expected 
to give rise to /nl of 0(1). This is clearly very close to what Planck will achieve, and it means 
that in order to be able to fully exploit the next new experimental results, the calculation of the 
non-Gaussianity induced by the non-linearities in general relativity and in the plasma physics needs 
to be done. 

Of course we are not the first ones to realize the importance of this issue. In particular, an effort of 
finding and solving the full set of second order equations is on-going by some groups |19| [20], [21"| [22] . 
For example, the three point function of temperature fluctuations in the particular limit where 
one of the wavelengths is much longer than the horizon at recomination and than the other two 
wavelengths has been consistently computed, finding an effect of order /nl ~ 1 (23] . This is still far 
from the required full calculation. The derivation of the second order equation requires to deal with 
some subtleties such as the formulation of the Boltzmann equation in curved spacetime that were 
neglected in former derivations such as the one of |19) . and that could affect the result at the order 
we care about. One of the two purposes of the present paper is to derive in a consistent way this 
set of equations. Unfortunately this will take us through a rather technical path, though we will try 
to make always clear the physical interpretation of what we are doing. We can already synthesize 
the main points of the derivation. We will work in the all-orders generalization of the Newtonian 
gauge, which is called Poisson Gauge. Writing Einstein equations in this gauge is tedious, but not 
particularly difficult. Writing down the Boltzmann equations in this frame is instead a little more 
subtle. The Boltzmann equation consists of a free-streeming term (so called Liouville term) which 
takes into account of the free evolution of the species one-particle distribution, and of a collision term, 
which takes into account how interactions among the particles affect the distribution. The collision 
term is expressed in terms of cross-sections that are measured and well expressible in Minkowski 
space. This means that, in order to write down the Boltzmann equation in Poisson gauge, one 
needs at each point to find the coordinate transformation to the local inertial frame, write down 
the Boltzmann equation there, and then transform back to the original frame. This will give us the 
recipe to write down the Boltzmann equations at second order in the perturbations. 

If schematically for the moment we write each perturbation 5 as 5 = 5^ + 5^ where 5^ is of 
order (5^) 2 , the resulting first and second order set of Einstein and Boltzmann equations will take 

1 This oversimplifies the structure of the possible non-Gaussianities [T3J [TJ, but this is enough for our 
introductory purposes. 
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the form of a system of equations of the form: 



L>[«5«] = , (2) 
D[8^} = S[5^ 2 } . 

Here D is a differential operator that is the same both for the first and the second order equations. 
The second order equation, unlike the first order one, has a source term S proportional to the square 
of the first order perturbations [^} 

Among the many terms that contribute to the source S, there are some (coming from the Comp- 
ton scattering collision term) in which one of the two first order perturbations is the perturbation 
in the number density n e of free electrons: 6 e = 5n e /n e . If one concentrates, as we will do, on the 
CMB, then 5 e affects the CMB anisotropy only at second order. This can be seen in the following 
way. The free electron density affects the CMB temperature through the Compton collision term, 
which determines the visibility function and the diffusion scale. The former gives the probability 
for a photon to originate from a given distance along the line of sight to the observer, while the 
latter accounts for the Silk damping at higher multipoles Vs. In a homogeneous universe, before, 
during and after recombination the radiation temperature decreases as a -1 , a being the scale fac- 
tor, irrespective of the electron density. Thus, a perturbation to the electron density changes the 
position of the last scaterring surface and the mean free path of the photons before decoupling, but 
not the observed radiation temperature. Therefore, since at first order we would perturb n e and 
keep the other quantities unperturbed, the CMB anisotropics are not affected by 5 e at this order. 
However, they clearly are so at second order. This explains why this quantity had never been com- 
puted before in detail. The first part of the paper will in fact be devoted to the computation of this 
quantity which is necessary to solve for the second order perturbations. As we will see, this task 
will not be completely straightforward because of the way recombination occurs in our universe. In 
particular, even in the homogeneous universe the collision term which controls the population of 
free electrons does not vanish (of course, otherwise there would be no recombination in the homoge- 
neous universe). This means that the same subtleties that we mentioned above about writing down 
the Boltzmann equation at second order will now apply also at first order. Further, even in the 
unperturbed case, recombination is not a very straightforward process, involving atomic transition 
between many different states. The treatment is simplified by the fact that matter and radiation are 
very close to equilibrium during recombination, with the most notable exception being the ground 
and first excited states of Hydrogen which are not mutually in equilibrium. This results in a series of 
approximations, which allow one to treat Hydrogen as an effective 3-level atom (ground state, first 
excited state and continuum) as was done in the classic work of Peebles [24J. In the end, one can 
write a single ordinary differential equation (ODE) governing the evolution of the electron density. 
Working in the same approximations, we shall derive the analogous equation for the perturbation 



of the electron density eq. (38). We also derive the equation for the perturbations of the matter 
temperature eq. ( p51~| 

2 The first order perturbations are non zero only once one assigns them some primordial non-zero initial 
conditions. 

3 One can try to directly perturb Peebles' ODE and write an equation for S e which then has to be comple- 
mented with an equation for the perturbations to the matter temperature, St m ■ The calculation of S e has been 
already attempted in |25j without accounting for the effects of metric perturbations and of the perturbation 
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At least on large scales, one can consider the primordial perturbation as representing a time 
delay St of the unperturbed solution in each point in space (of the order of W~ 5 H^ 1 ). In this case, 
5n e ~ h e 6t, and since the time-scale of recombination is rather quick, h e /n e ~ 15 H (that is 5 times 
the expansion rate), we expect 5 e ~ 5 x 10 -5 , an enhanced perturbation with respect to the standard 
ones. This is in fact what will be the result of a detailed calculation. 

Because of this enhancement, it is natural to wonder if this might affect the CMB, for example 
by generating an /nl of order 5, which could be potentially detectable by Planck. We find the way 
S e changes the CMB is mainly by delaying the time of recombination, by modifying the shape of 
the visibility function, and by changing the diffusion scale. The result is that 5 e gives a potentially 
detectable effect equivalent to /nl ~ 5. Since the calculation and the interpretation of this result is 
not completely straightforward, it will be the subject of a companion paper [28 . 

The paper is organized as follows. In sec.[2]we review recombination in the unperturbed universe. 
In sec. [3] we derive the first order equations for recombination in the perturbed universe and discuss 
the results of their numerical integration. In sec. [4] we derive the full set of second order equations for 
cosmological perturbations. In sec. [5] we summarize our results. In App. [A] we review the formulation 
of kinetic theory in curved spacetime. In App. [B] we derive in a rigorous way the perturbations to 
the escape probability. In App. [C] we give an erratum to the Compton collision term at second order 
which is already present in the literature. 

2 Review of standard recombination 

Here we review the standard Peebles' [21] treatment of recombination in the homogeneous universe. 
We will also state what are the approximations that are implicitly done in this treatment. When 
dealing with cosmological perturbations, there will be additional timescales corresponding to the 
frequency of the various modes, and we will have to check in which regime the standard approxima- 
tions still hold. For simplicity of discussion, we will neglect Helium, since its treatment is completely 
analogous to the treatment of Hydrogen recombination. We work with units in which the speed of 
light, c, is set to unity. 

2.1 The effective 3- level atom 

Recombination in the early universe is an example of the so called Case B recombination, i.e. all 
levels with n > 2 are in radiative equilibrium, but the Lyman a line is optically thick. This 
means that a 2P — > IS transition generates a high energy photon which immediately causes a 
IS — > 2P transition in a nearby atom, which results in an overpopulation of the n = 2 level 
relative to radiative equilibrium. The recombination to the IS level proceeds in two channels: either 
through the redshift of Lya photons out of the absorption line, or through a two-photon decay 
from the 2S level. Thus, in the homogeneous universe one can treat Hydrogen recombination using 
the effective 3-level atom picture described by [21]. Due to its simplicity, it is tempting to use 

to the escape probability, and in [26] |2Z] , where only a non rigorous derivation of the perturbations to escape 
probability was given. We will show that the formulas given in [55] [57] are indeed correct apart for irrelevant 
corrections. 
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the Peebles equation for recombination in the presence of perturbations. However, the Peebles 3- 
level atom is self-consistent only because the physical timescales involved in recombination span 
a huge range of scales. Thus, a number of processes can be regarded as instantaneous allowing 
one to make several quasi-equilibrium approximations which we will discuss in detail below. In the 
presence of perturbations, however, one should consider whether the perturbation timescale is long 
enough to allow for the same approximations used in the homogeneous case. In order to answer this 
question, we will list the most important approximations together with their relevant timescales, 
where appropriate, which allow one to calculate the recombination history to an accuracy of about 
1% |29j . We will delay the discussion of the approximations allowing for the calculationg of the 
Lya escape probability to the next section. As we go along, we will gather all relevant timescales in 
Figures [l] and [2] For reference, the recombination timescale is t rec = \n e /h e \ ~ 5 x 10 4 yr (= 15 Mpc 
in conformal time) during recombination (cf. Fig. [TJ. 

The first relevant timescale enters in the calculation of the photoionization and recombination 
cross sections. This calculation relies on the fact that the electrons and baryons obey the Maxwell- 
Boltzmann distribution with a common kinetic temperature, Tm- This is justified beacuse the 
thermalization timescale is given by the collision timescale, which is always less than about a year 
for the relevant redshifts (z > 500), corresponding to a comoving k > lOkpc -1 (see Fig. [l| which 
is way beyond the scales we care about. Other processes, such as collision, stimulated emission and 
stimulated recombination could be in principle important for recombination, but they are not, as 
shown in [29 . 

Next, another important assumption is that all levels with n > 2 are in radiative equilibrium. 
This comes from the fact that, apart from the Lyman series, the number of photons generated by 
atomic transitions can be completely neglected with respect to the number of photons in the Planck 
spectrum at the corresponding frequencies [30 a . Since the recombination rates involve a collision with 
an electron, they are orders of magnitude less frequent than radiative transitions, which therefore 
should keep the populations with n > 2 in thermal equilibrium with the radiation. Thus, the ratio 
of the populations of two excited states i,j>2 is given by the Boltzmann factor with a temperature 
given by the radiation temperature Tr entering in the black body distribution: 

"1 = ii e iEi- E i)/kBT R ^ / 3 n 

where gi and Ei are the degeneracy and ionization energy of the i-th level, respectively, and ks is 
Boltzmann constant. This assumption has been checked in [29], where it was shown that the above 
approximation is valid for levels which are nearby in energy (e.g. 7170/7175), but for z < 1000 it fails 
to reproduce the ratio of populations for levels which are well-separated in energy (e.g. 7170/712) 
due to the photoionization becoming inefficient. Nevertheless, this induces an error only at lower 
redshifts, z < 800, which amounts to no more than about 10% of the electron fraction x e since at 
these redshifts the populations of the highly excited states are very small [29 . This correction can 
be accounted for in the Peebles 3-level atom by the introduction of a fudge factor to the absorption 
rate |29j as discussed below. However, this tells us that the timescale for relaxation among very 
different states is not much faster than the recombination time-scale at z ~ 800. Fortunately, the 
photoionization rate grows very rapidly as we increase the redshift (see Fig. Q), and therefore, even 
though we are interested in perturbations whose timescale is faster than recombination, we expect 
the thermal equilibrium approximation among the states n > 2 is a good approximation to ~ 10% 
for momenta k up to order 1 Mpc -1 . 
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Figure 1: We summarize all important timescales related to recombination in this figure and 
m Fig. |J Here we concentrate on the timescales that control the 3-level atom approximation. 
The peak of the CMB visiblity function is at z = 1090. The expansion timescale is defined 
as (37i) _1 ; the Lya relaxation timescale is given by eq. (29) of [H]; the collisions timescale 
corresponds to atomic collisions, and is irrelevant at early times; the Kompaneets timescale 
corresponds to the timescale of Compton heating. Of the timescales given in this figure, 
the Lya relaxation imposes an upper limit of k ~ 1 Mpc -1 on the validity of the perturbed 
Peebles equation. 
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Given the above approximations, we can now proceed along the lines of [24], and write down the 
equations for the evolution of the electron number density n e and the population, n 2 , of the first 
excited state: 

+ 3Fn e = -AR C _ 2 , (4) 

^ + 3Hn 2 = AR C ^ 2 - AR 2S ^ 1S - AR 2P _> 1S , (5) 

where H = d\aa/dt is the Hubble parameter and a is the scale factor; AR 2 s->is is the net rate 
of 2S — * IS two-photon transitions, AR 2 p^,i$ is the net rate of IP — * IS transitions, and Ai? c _>2 
is the net transition rate between the continuum and the first excited state, including transitions 
which involve intermediate states such as continuum — ► (n > 2) — ► (n = 2). Notice that we did not 
include transitions from n > 2 to 15 because they are exponentially suppressed. Because the upper 
levels are in radiative equilibrium, one can write AR C ^ 2 in terms of a total recombination olb and 
photoionization (3b coefficients as follows 



AR C ^ 2 = a B (T M )n e n H u - (3 B {T R )n 2S , (6) 

where nnn is the number density of free protons, which after neglecting Helium, reduces to n e = 
nnn- Neglecting stimulated recombination, as can only depend on the matter temperature (and 
negligibly on the number density of electrons [3T] ) , since recombination is a collisional process 
which depends on the kinetic temperature of the electrons Tm • The Case B total recombination 
coefficient is given in [31], and an adequate fitting formula for oib(Tm) is given in [32]: olb(x) = 
1.14 x 10~ 13 aj/V(l + cy d ) cm 3 s , where y = x/10 4 K, x is the temperature in Kelvin degrees, 
and a,b,c,d are dimensionless constants of order 1 (a = 4.309, b = —0.6166, c = 0.6703, d = 
0.5300), and 1.14 is the fudge factor introduced by [29] to account for non-equilibrium populations 
of the excited states. The total recombination coefficient is given by a sum over the recombination 
coefficients, aj(T/vf ), to each excited state. One can derive a total photoionization coefficient (3b from 
detailed balance considerations (see e.g. [29 J. In local thermodynamical equilibrium (LTE), using 
equilibrium populations denoted with a superscript eq one can write Pin^ q = n^n^^ai, where PiiTg) 
is the photoionization coefficient for each excited state. Then, the total photoionization coefficient 
from the first excited state is given by /3b = Y2i>i A(7fi)(^i / 'n 2 ) eq , which can be combined with the 
LTE expression and ([3| to give ([7| below. Since both the individual photoionization rates and the 
population ratios ^ depend only on Tr, the total recombination coefficient can also depend only on 
Tr. For the purposes of comparison with the Peebles 3-level atom, [29] obtains an expression for (3b 
which depends only on Tm- During recombination, this is an adequate approximation since Compton 
scattering keeps Tm equal to Tr, even for the faster timescales associated with the perturbations. 
We explicitly verify this by keeping Tm and Tr distinguished and by solving for both of them. In 
this case, (3b(Tr) is given by: 

Pb(Tb) = a B (T R )e-^« i^^ 12 > (?) 

where hp is Planck's constant and m e is the electron mass. 
The 25 — > IS net transition rate is given by 

Ai?25^i5 = Ms^isn 2 s - A 1S -> 2 s{TR)nis , (8) 
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where A 2 s-*is = 8.22458s 1 |33j. From detailed balance n^ s Ais-^2S = ^s^-ss^is, we obtain 
Ais^2s(Tr) = A 2S ^i S e- B ^ kBT K, where B 12 = E 1 - E 2 . 

We will derive the net rate for the 2P — ► 15 transition (28) in the next section. Neglecting 
stimulated emission, and using the fact that the Lyman a line is extremely optically thick we obtain 



n 2P -3n ls e- Bl2/kBTR ] , (9) 



where A 2 \ is the rate of spontaneous emission from 2P — ► 15, and P = P 2 \ is the escape probability 
for the Lya photons discussed in the next section. Notice that, if for a moment we set P to one, 
the above expression is the standard rate 2P — ► 15. As it will become clear in the next section, 
the multiplication by the probability P takes into account of the fact that just a fraction of these 
transitions are effective. The rates R 2 p->is /"-hii and Ris~, 2 p /rami per HII are shown in Fig. [2j We 
can see that Lyman a escape is important only very early on (z > 1400), and that the two-photon 
decay dominates most of the recombination history. 

We can now begin to solve the system of equations Q and ([5]). First, we can neglect the left 
hand side of §5§ since the recombination and Hubble timescales are large compared to the net atomic 
transitions timescale per excited atom which is about A^o^ig ~ 0.1 s for each of the three net rates 
appearing on the right hand side of ^ . This is equivalent to saying that the n 2 s and n 2 p reservoirs 
are in quasi-equilibrium [33]. In a perturbed universe, this approximation is again justified, since 
it corresponds to neglecting the perturbation timescale with respect to the net atomic transitions 
timescale. Assuming equilibrium populations of n 2 p and n 2 s due to fast radiative transitions gives 
n 2 = 4n 2 p/3 = 4n 2 s- Thus, combining (|5]),([6]),([7j),(|8j an d ([gj^ we obtain an algebraic equation 
which can be solved for n 2 , obtaining: 

gl2 

a B n 2 e + (A 2S ^is + 3PA 21 ) e fe B T «(n b -n e ) 
n2=4 o Pi T a TTi ' ( 10 ) 

where is the number density of protons (both in ionized and atomic Hydrogen). Plugging this 
solution for n 2 into Q, we obtain the standard Peebles 3- level result, with the distinction between 
Tp and Tm made clear: 

din 

+3Hn e = Q[n e ,n b ,T M ,T R ,H] , (11) 

where 

Q = -[a B (T M )n 2 e - /3 B (T R ){n b - n e )eT B ^ k ^\C H , (12) 

and 

r =1 Pb{Tr) 

H 3PA 21 + A 2S ^ 1s + /3b(T r ) ■ { ' 



As we will show in the next section, 



3 PA 21 « , (14) 

rib - n e 



where v 2 % is the frequency of the Lya photons. Here and in eq. (10) we have approximated nxs ~ 
rib — rie, since even in the presence of the Lya bottleneck, the ground state population dominates the 
rest of the excited states. Note that in the above derivation we could not neglect the left hand side of 
Q as we did in (O, since the net atomic transitions timescale per electron is A^g^g X n e /n 2 ~ t 
(see Fig. [2J. Thus, in the presence of perturbations, the left hand side of Q will receive non- 
negligible corrections, as we will show in Section [3] 



rec 



s 




Figure 2: The comoving rates per HII ion governing Hydrogen recombination. The recom- 
bination rate is defined as \n e /n e \. Here we concentrate on the rates that affect directly the 
recombination rate. 



9 



2.2 Treatment of Lya photons 

In the expanding universe photons are "advected" in momentum space by the cosmological redshift. 
The photons blueward of the Lya line are redshifted into the line, where they diffuse until they reach 
the red wing of the line from which they can "escape". In this section we derive the Boltzmann 
equation for the evolution of the photon distribution function. Solving that equation will allow us 
to find the escape probability of Lya photons and derive eq. ^ . 

In general, the physical photon number density in each frequency bin of width Sv is given by 
8nu 2 F-y(i'(t), t)5u, where 8irv 2 gives the number of states per unit volume per unit frequency range, 
and F-y(u,t) is the photon phase-space distribution (one can think of it as the occupation number 
of each state). Note that we have included the photon degeneracy in 8iru 2 . Following a given 
frequency bin as it is being redshifted, the number of states in a comoving volume V is given by 
a? (t)8itV v{t) 2 5v(t) , which is constant in time. Therefore, the rate at which the number density of 
photons increases due to photons which are redshifted into a bin of a given fixed frequency v between 
a time t — St and t is simply given by 

R in = 8itv 2 5vF 1 [v{1 + H5t),t- St]/ St , 

where we used Sv/St = —Hi/. The rate at which photons are redshifted out of the bin is instead 
given by 

Rout = 8 , kv 2 SvF~ 1 \v, t\/5t . 
The net flux divergence equals therefore 

Rin- Rout = -8w 2 (^-Hv^)sv = -8w 2Cl ^Sv. (15) 



dt dv ) dt 

The above expression should equal the net rate at which photons are being absorbed in a certain 
line representing a transition from the state j to the state i: 

net absorption rate = —[AjiUj — (Bijm — Bjinj)2u 3 F^{y)\cf){ii)5v , (16) 

where Aji,Bji and Bij are the standard Einstein coefficients, and §(y) is the line profile In the 
above equation we used that the intensity of the radiation field is given by I = 2z^ 3 i ? 7 , and that 



In the photon collision term in the above eq. ( 16 1 we assumed that complete redistribution holds by 



setting the emission and absorption profiles equal. In reality this is far from true in general processes. To see 
this, consider a 2-level atom with levels 1 and 2, and for simplicity let us concentrate on the Lya transition 
(i.e. level 1 is sharp) which is the most important during recombination. When an atom is excited to level 
2, it spontaneously emits a photon within ~ 1CP 9 s. During this whole process we can assume that the atom 
moves with constant velocity, since collisions and atom recoil can be neglected. In its rest frame, the atom 
emits a photon with the standard Lorentzian profile. This is because we can assume that the absorbed and 
emitted photons are not correlated in the atom rest frame since the absorption and emission events are well 
separated in time. However, the atoms follow a Maxwell-Boltzmann distribution, which implies that in the 
laboratory rest frame, the absorbed and emitted photons are correlated, with the correlation given by the 
so-called redistribution function R{v' ', h'\ v, h) which depends on the frequencies of the incoming and outgoing 
photons and on their directions. Neglecting angular redistribution |35j . this results in an emission profile given 
by 4>2i{y) = J dv 1 1(y')R\\{y' , v)j J dv'I(y')(j)i2{y'), where R\i{v' ,v) is the partial frequency redistribution 
function subject to the above set of assumptions |36j . To understand how one obtains the redistribution 
function for an n-level atom, one should start from the formulation of the radiative transfer problem in terms 
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collisions are unimportant during recombination as shown by [29 . Equating the net absorption rate 
to the flux divergence we obtain the following differential equation 



dFy 
'~d~t 



1 



iirv 



[A 



3l n,j 



2v i {Bijni - Bjirij)Ey]<l>(i/) 



(17) 



This is exactly the Boltzmann equation for the photon distribution. In App. [A] we review the 
covariant formulation of the Boltzmann equation necessary for perturbing it. 

Before actually solving this equation, we can already estimate the rate Ait^p^is, which must 
equal R ou t — Rin- We can approximate eq. (15) by AR 2 p^\s ~ 8ttHv 3 [F^(vji) — F 7 (i>g)] with vr 



and vb two frequencies just to the red and blue of the line. Just blueward of the line, the spectrum 
converges to the black body spectrum B V: i.e. 2v b F~ 1 {vb) = B v . For optically thick lines, the photon 
distribution function saturates to S/2V 3, at vr, with S being the standard radiative transfer source 
function (see (20) below). Thus, we can write AR2p~,is ~ 4:7rH(S — B u ), which as we will show 



below (cf. equations (26) and (28)) is an excellent approximation for optically thick lines. By using 
the standard relation Aji = g^v 3 Bji, it is easy to see that this result is equivalent to eq. (JoJ) that we 
used in Section 12.11 

Let us now proceed and solve the Boltzmann equation (|17|) to obtain the rate AR C ^2- Following 



we can reparametrize (17) using the proper frequency v(t) of an individual photon as a time 



parameter, since v(t) oc a(t) is monotonous: 



7" 



dFry 

dv 



where 



S 

2^3 



4irH 



F-, 



<f>(u) 
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BijTli Bj {Tlj 



is usually referred to as the Sobolev parameter, and 

S = 



Ajijij 



BijTli BjiTlj 



is usually called the source function. The solution to (18) is given by 



F y (u) = F^v B )e- T ^ + e 



Jug 1V 16 dv' 



with 



r{v) 



dv' 



(18) 



(19) 



(20) 



(21) 



(22) 



of quantum mechanical generalized redistribution functions. For a systematic discussion of how this is done, 
see |37j and [38] - Using the correct redistribution function Ru for the Lya transition, the radiative transfer 
equation can be parametrized in the form given in [35 for example. It was shown in [35] that the line profile 
averaged intensity, which is what enters the net transition rate, is weakly dependent on the choice of the 
redistribution function for 7 (defined below in eq. (19)) in the range 10~ 4 to 1. Outside this range, the result 

i(v)4>(v')) is practically identical to the one using the 



using complete redistribution (corresponding to i?cR ; 
Ru redistribution function. For the Lya photons 7 is of order 



10 /(l — x e ) during recombination. Thus, 



once we have non-negligible recombination, we can work with the complete redistribution approximation. 
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The optical depth t{v) is calculated along incoming rays. We will see later that the final result is 
independent of whether one considers incoming or outgoing rays. The choice of vb is arbitrary, but 
for the approximations below to work, we choose it to be just blueward of the line. The frequencies 
just blueward of the line are not affected by the line, and therefore, we can write Fj(vb) = F^vb) ~ 
Fbb{ v ij)i where is the Bose-Einstein (black body spectrum divided by 2v 3 ) distribution. 
The net transition rate from j — > i is given by 

oo 

ARji = J dvcj){v)[Ajinj - 2u 3 (Bijni - Bjin^F^] . (23) 
o 

The integrand is multiplied by the line profile which is centered around Vij and has a width Av given 
by the thermal broadening of the line. Since v plays the role of a time coordinate, in configuration 
space this implies that the contribution of the number densities entering the integrand will be only 
from a spacetime patch of scale L$ = v thermal/-^ which we call the Sobolev scale [^J Here v thermal is 
the average thermal speed of the atoms. We can take the Sobolev scale as approximately constant 
during recombination and equal to L$ = 6.3 kpc (~ 20 yr in physical time). This is very small 
compared to the recombination timescale (which is also the timescale of variation of the source 
function), about 15Mpc. 

This allows us to make the so called large velocity gradient or Sobolev approximation jlQ] where 
we neglect the spatial and temporal variations of the number densities over the Sobolev scale. This 
is not an innocuous approximation: this amounts to neglecting the relaxation time over which F^ 
becomes quasi-static in the line. The validity of this approximation has been checked by |41| . where 
they show that the relaxation timescale of Fj, defined in |TT] as the timescale over which the intensity 



in the line reaches within 10% of the quasi-static value, is t re n ~ 10 3,5 yr in physical time (t T ^ ~ 1 



comoving Mpc ) at the peak of the visiblity function (see Fig. |lj). This effect, which becomes 
relevant for modes of order IMpc" 1 , is not very important because the transition 2P — > IS is also 
not very important 



Neglecting variations of the source function and the Sobolev parameter we can rewrite ( 23 ) as 



oo 

ARji = AjiUj - (Bijm - BjiUjpvfj J F^{y)dv . (24) 



o 



In the same approximation we can take 7 outside of the integral in (22), and obtain t(u) = 
7~ 1 (fj ? ') f^ B dv'(j)(v'). At this point we are allowed to extend the integral range to vb 00, since 
the contribution of the line profile at these frequencies is negligible. Thus, we can finally write 



in the Sobolev approximation. The integral in (24) can be evaluated using (21) by a change of 



For a single photon, the Sobolev patch extends a scale L$ in the time direction, and a much smaller scale 
in the space direction, since the photon experiences a random walk. The relevant scale is given by Lg, the 
largest of the two. 

6 The effect of the Lya escape probability is about 10% for k ~ 0.1 Mpc" 1 (equivalent to a CMB multipole 
I - 1500) and about 25% for k ~ 0.2 Mpc" 1 (/ ~ 3000). 
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variables using dr 



[y)dv. In the end we obtain 



oc 

2v% J F^{v)dv = P 3i 2vlF bb {vij) + (1 - P Si )S(vij) , 



(25) 



where 



P, 



j 1 



dv(j){y) exp[— t{v)\ = 7 ^1 — e l /<\ 



(26) 



The prefactor Pji is called the escape probability, as it is the probability for a line photon to escape 
to infinity without scattering, once it is spontaneously emitted. To see why this is so, let us consider 



the different components of the line-averaged intensity given by the above equation ( 25 ) . The first 



term, 2vf-PijF bb {vij) : is called the direct component, since this is the intensity coming from the 
blue part of the spectrum which equals 2uf-F bb {vij) times the probability for it to "escape" to the 
center of the Sobolev patch. The second term is called the diffuse component of the radiation field, 
(1 — Pji)S(vij), since it is that part of the intensity which is locally produced, and which does not 
escape the Sobolev patch. Using the definition of the source function, we can write 



~\~ fljSjj](l Pji)S — TljAjiP, 



(27) 



which means that the net rate of producing local (diffuse) photons equals the rate of spontaneously 



emitting escaping photons. Using (|20j), (24) and (26) we can finally write 



ARji = PjilAjiUj - {Bijm - B ji nj)2vf j F bb (i'ij)] . 



(28) 



Neglecting spontaneous emission, and specializing to the Lya transition, we recover ([9]) as wished. 
We notice that in the Sobolev approximation, after averaging over the line profile (i.e. integrating 



by J du(f>(u)), eq. (21) becomes: 



F 7 = PF bb + (1 - P) 



1.1 J 



(29) 



where the bar stays for frequency averaging. We will find a very analogous expression in App. [B] 
when we will generalize to the perturbed universe. 

As anticipated, we also notice that the Sobolev approximation ensures that the forward and 
backward escape probabilities are equal, since 



P» 



dv(f){y) exp 



-7 Vu) 



dv c 



du<p(y) exp 



-7 Vij) / di/(j)[y') 
Jo 



In the discussion above we neglected the secondary distortions to the radiation spectrum [29] . 
These distortions come from the redshifting of line photons into a neighbouring line. This effect is 
trying to invalidate the approximation F^(vb) — Fbb{ v B) that we used in the derivation of the escape 
probability. This approximation holds everywhere outside the Lyman series, where the relative 
intensity distortion is AI j '{2v^ F bb ) < 10 -6 [30]. The secondary distortions due to the higher-energy 
Lyman lines on the lower-energy Lyman lines has been assessed in [32]. The biggest effect comes 
from the secondary distortion caused by the Ly/3 line on the Lya line. However, the number of extra 
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Lya photons per H atom due to this effect is 0.2% |42j . therefore the effect of perturbations on the 
number of redshifted Ly/3 photons will have negligible impact on 8 e . 

This concludes our discussion of Hydrogen recombination in the unperturbed universe and of 
the escape probability. More details on the Sobolev escape probability can be found [35], while a 
recent review on of recombination is given in [13] and the references therein. 

Throughout this discussion so far we have neglected Helium recombination. We can easily 
reintroduce it: in this case, the Peebles equation for n e becomes an equation for nnn, and one ends 
up with two extra equations for the abundance of Hell and Helll, given in a convenient form for 
example in [25J R 



3 Perturbing recombination 

We are now ready to derive the equations for 6 e and STm- Throughout the rest of this paper, we will 
work in a spatially flat FRW background. We work with a metric with positive signature. Greek 
letters denote spacetime indices; latin letters denote spatial indices. In this section, we will work to 
first order in the metric perturbations in the Newtonian gauge in which the line element is given by 

ds 2 = a 2 (r/)[-(l + 2V)dr) 2 + (1 - 2$>)dx 2 ] , (30) 

where rj is conformal time, and <1? are the Newtonian gauge potentials, and a(rj) is the scale factor. 
A dot will denote a partial derivative with respect to ij. 



3.1 Electron density perturbations 

In the previous section we found that all the approximations entering in the Peebles 3-level atom 



equation (11 ) are still valid in the presence of perturbations with k < 1 Mpc . This limit comes both 
from the timescale over which the excited states reach thermal equilibrium, and from the timescale 
over which the photon distribution reaches quasi-equilibrium in the Lya line. Limiting ourselves to 
k's lower than this, we are therefore going to study the Peebles 3-level atom in a perturbed universe. 
In order to do this, we have to derive the escape probability in the presence of metric perturbations, 
and we also have to write down what the Boltzmann equation for the electrons is in curved spacetime. 
The last step is really necessary because, contrary to the treatment of first order perturbations, in the 
case of recombination the collision term does not vanish in the homogeneous universe. This means 
that we have to consistently treat the metric perturbations both in the Liouville (or free-streaming) 
term, and in the collision term. 

We start from the covariant form of the Boltzmann equation. The covariant formulation of the 
Boltzmann equation in a curved spacetime is reviewed in App. |Aj to which we will refer on several 
occasions in order to use some results from covariant kinetic theory. From (110) in App. |A| the 
Boltzmann equation can be written as 



7 Note that we do not need to revert to the more detailed equations given in [43 , because Helium recom- 
bination has little effect on the CMB. 
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where f(x l ,Pj;r]) is the one-particle distribution function; P^ = dx^/dX where A is the affine 
parameter chosen so that P^ is the momentum of the particle (dX coincides with the time interval of 
a local inertial observer divided by the energy observed in the same frame: dX = dti nert iai/-Einertiai); 
C[f] is the collision term given in (114 ) in App.[A|. The one-particle distribution function f(x^, P v , 77) 
is a scalar function of the phase-space variables, which gives the distribution of particle momenta at 
every spacetime point. A very simple way to understand the structure of the Boltzmann equation 



is to notice that the left hand side of eq. (31 ) is just the total derivative of the distribution function 



with respect to the affine parameter A. Since these are scalar quantities, the collison term C[f] on 
the right hand side is the same as the one measured in a local inertial. The factor of P° is due to the 
gravitational redshift in converting the coordinate time to the time of the local inertial frame where 
C[f] gets a very simple form (see App. |A| . In fact, it is very useful to express the collision term at 
a certain given coordinate x^ 1 using the momenta p 1 of the particles measured in the Local Inertial 
Frame Instantaneously at Rest with respect to the Comoving Observer (LIFIRCO) at fixed x l . In 
this case the cross sections have the same expressions as in Minkowski space, and the collision term 
contains no metric fluctuations. All the metric fluctuations are confined to the Liouville term. This 
leads us to often use mixed coordinates (x^p 7 ), which we call "nice" coordinates, where the spacetime 
coordinates are generic, for example they can be the ones associated to the Newtonian gauge, while 
momenta are the ones corresponding to the LIFIRCO. The distribution function expressed with 
these nice coordinates we denote with F(x l ,p J ;r]) = f(x l , Pfc(x*,p 7 ; 77); 77). 



From the Boltzmann equation (31), one can derive the conservation of 4- momentum for a given 



particle species given in JTT3] ) in App. |A} 



(n U^) ; 



C[f] , 



(32) 



where n is the number density of particles measured in a local inertial frame momentarily at rest 
with respect to the fluid, tt is the invariant measure in momentum space on the mass-shell, which 

p° is the particle 



reduces to g& eg d?p/E in a local inertial frame (see eq. (102) in App. [A|, where E 
energy as measured in that frame, and g deg is the number of internal degrees of freedom for each 
particle of the given species. We have defined the 4- velocity in such a way that the particle 
current can be expressed in the following form 



1 1 



We can write equation (32) in nice coordinates: 



d 3 p 
~E 



C[F] . 



(33) 



The Newtonian gauge 4- momentum P^ is related to the 4-momentum p^ = [E,p l ), measured in the 
LIFIRCO, via the tetrad (103) discussed in App. [Aj We call the 3-velocity of the fluid measured in 
the LIFIRCO as v 1 = p l /E. At first order, the relation between U^ 1 and v % can be found in e.g. |44j . 



or by using (90), which we will derive later. To first order this is given by: 

U° = , 



[1 - #) , u l 



(34) 



For a nuclei species n, number conservation holds (i.e. the integrated collision term in (33) 
vanishes in this case), since nucleosynthesis is long over by the time of recombination. Therefore, 
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from ( 33 ) and ( 34 ) we have 



7=4°) + 3n^Ti = , 
9 n + S n - 3$ = , 



(35) 
(36) 



where we have expanded n n = rin (1 + S n ), and we have defined 9 n = v l ni (this coincides at first 
order to the 9 variable in [2]). The above equations are the usual continuity equations written in 
the Newtonian gauge (see e.g. 



For the electrons the collision term does not vanish, and we obtain, from (33): 

n(°)+3ftn( 0) =«Q (0) , 

nf\5 e - 8 b ) = a(*QW + 5Q) - 6 e aQ<® , 



(37) 
(38) 



where TL = a/ a = a H . In deriving (38) we used the baryon number conservation equation (36) to 



express in terms of the perturbation to the baryon number density 5 b . We have split the integrated 
recombination collision term Q into a zeroth and a first order part as 



Q 



nf\l + 5 e ),nf\l + 5 b ),T$ + ST M , T R 0) + 5T R , + S H ) 



Q {0) + 5Q . 



(39) 



Notice that, as we have explained before and we explain more in detail in App. the collision term, 
being a local and scalar quantity, once expressed in terms of nice momenta p l, s does not contain 



metric perturbations. This means that it is the same function as in the homogeneous case (eq. (11 )) 



expressed in terms of the perturbed quantities rie (1 + S, 



n 



(0)/ 



(1 + 5 b ), Tj? + 5T M , and T£> + ST R . 



(0) 



The only exception to this reasoning is the escape probability, which, as we saw in the former section, 
is rather a non local quantity. However, for perturbations with wavelength much longer than the 
Sobolev length, we can imagine to recover the expression for the perturbed escape probability from 



the result in the homogeneous case. In fact, by inspection of eq. (14), we see that what controls the 



escape probability is the Hubble rate if, which in a homogeneous universe is exactly equal to 
the local baryon velocity divergence, [7^ /3 = H. This is not a coincidence, as in fact what 
controls the probability for a photon to escape is the redshift of the neighboring atoms, which 
comes from their relative velocity. It is at this point easy to guess what is the expression for 
the escape probability (and therefore for the integrated collision term Q above), in the case 
of a perturbed universe: in the expression for Q, we need just to replace H with £7^ /3, or 



equivalently, as already expressed in eq. (39), H with H(l + 5h), where 
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H 



3H 



m 



(40) 



By using these arguments, the above expression for the perturbed escape probability had 
been derived in App. D in [27] . A rigorous derivation of the perturbations to the escape 
probability is given in our App. [B]by angle averaging the photon distribution function using 
the Boltzmann equation. We show that there are extra terms that appear in the expression 
for 5h due to the change of the baryon fluid velocity due to the acceleration by pressure 
forces and Thomson drag over a time interval ~ Lg, but we argue that those contributions 



are negligible. This gives a rigorous, though more convoluted, derivation of eq. (40) 
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So far we have neglected Helium recombination, which can be introduced by complete 
analogy. The perturbations to the Helll density are given by perturbing the Saha equation, 
since Helll recombination procedes very nearly perfect local thermal equilibrium. The equa- 
tions for the HII and Hell fractional perturbations are given by (38) with e replaced by HII 
and Hell, respectively, and Q replaced by their respective collision terms which one can read 
off from Section 3 of [35] . 



3.2 Kinetic matter temperature perturbations 

In this subsection we derive the equation for the perturbation of the kinetic matter temper- 
ature. By matter, we refer to the collection of electrons and H and He ions and atoms in all 
energy levels. Collisions and Coulomb scattering are very fast, and force all these components 
to have the same temperature and the same velocity 

The matter temperature measures the thermal kinetic energy of the matter species, there- 
fore, we look for an equation which includes all important effects which can change the kinetic 
energy of the matter species. Apart from the usual momentum redshifting, those include elec- 
tron Compton scattering, bremsstrahlung, and other subleading interactions. Both Compton 
scattering and bremsstrahlung keep Tm and Tr in thermal contact, but since bremsstrahlung 
has a timescale about 10 2 longer than Compton scattering, we can safely neglect the free-free 
transitions. 

Another effect comes from the fact that during recombination the actual number of free 
particles changes. In thermal equilibrium, each particle carries on average an energy equal 
to 3/c^T/v//2. When an electron combines with a proton, its kinetic energy is lost to radia- 
tion, and now there are fewer particles to share the remaining energy. This effect is usually 
referred to as molecular weight change, and it is entangled with another effect, usually called 
recombination-cooling/phoionization-heating, which takes into account that the kinetic en- 
ergy of the electron that recombines might actually be higher or lower than the average, which 
means that the remaining particles will have respectively a lower or higher average kinetic 
energy, and therefore a lower or higher matter temperature. These effects have been shown in 
[29] to have a negligible effect in the unperturbed universe. This is true even in the perturbed 
case, though we keep the molecular weight because the derivation becomes more transparent. 
Finally, the bound-state — > bound-state (bb) transitions affect the kinetic energy just through 
a negligible recoil energy, and we can therefore neglect them. 

We can find an equation for the matter temperature by defining a kinetic stress energy 
tensor for all the matter, Tm^v and writing down its divergence. The word kinetic stays 
simply to mean that we are neglecting the differences in potential energy (as we said, this is 
ok because the bb transitions do not affect the matter kinetic energy). Being a perfect fluid 
to a very good approximation, we can write the kinetic stress energy tensor in the standard 
form 

= (pm + Pm)U»U" + <r PAf , (41) 
where pu and pm are the energy density (without takeing into account of the electron potential 
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energy) and pressure of matter measured by a local inertial observer instantaneously at rest 
with respect to the fluid. They are given by (setting ks = 1 from now on): 

Pm = [nf(m p + m e ) + n%\{m He + 2m e )](l + 5 b ) + -p M , (42) 
(o) 



p M = n t (T ( M > + 5T M ) 



with 



ni 0, =ni°)+np), 

n t = n£\l + 5 b )+nW(l + S e ) . (43) 

Here n p counts all protons and H atoms (in all excitations), counts all He ions and atoms. 
In the above we used that the H and He have one and the same overdensity, 5 b , which follows 



from the first order continuity equation (36) and the fact that they have a common bulk 
particle velocity. 

In order to write the divergence of this stress-energy tensor, we use the fact that for each 
of the species 

d 3 P i nUT1 f d 3 p 



T 1 = Sde g / -rprV^P^Puf = 2dc g / -fP'PuF , (44) 



where in the last passage we have used eq. (104) of App. [Aj We can write the covariant 



derivative of T^ v of any species using the distribution function for that species [46J 

where A is the affine parameter defined as in the former subsection. The above equation is a 
differential geometry identity jl6] for the tensor given by (44), i.e. it holds for any species. 
By using the geodesic equation DP^/dX = P u P' 1 l/ = 0, we can express this in terms of the 
collision term 

Tr v = J 7T = J 7T P»C[f] . (46) 

We can apply the above identity also to the sum of all the matter components. In this 
case, on the right hand side there would be an integral of the electron Compton scattering 
collision term (the proton Compton scattering is negligible) , and of the collision term due to 
bound-bound and bound-free transitions. As we said, though associated to the emission of 
even highly energetic photons, these last ones only affect the evolution of the population of 
the excited states, or equivalently of the potential energy of the bound electrons as it gets 
diluted into radiation, but they do not affect the matter kinetic energy^} Therefore, for what 



3 The recoil energy due to a recombination event is very small. 
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concerns the kinetic part of the stress energy tensor defined above, its divergency is given by 
simply considering the Compton collision term. 



rjiflU 

1 M \v 



(47) 



After neglecting polarization, the second order Compton collision term is given by [UJ 
by expanding the Klein-Nishina cross section to second order in the baryon velocity (first 
order in Tm)- Their expression is further reduced to a useful form in [IH] (BMR1 from now 
on), their eq. (4.42). In App. [Cj we give an erratum to the second order collision term given 
in both papers. Using the corrected expression for the Compton collision term, the \x = 
component of the above integral can be evaluated transforming the momenta to the ones of 
the local inertial frame at rest with the rest frame of the baryon fluid: 



e ,y 



1)5=0 frame 



C 7e [F 7 ]d 3 p 7 



1 - 



-A 



c 



(48) 



where 



Ac = n e T R (T R - T Mj 



(49) 



and e 11 ^ is just the tedrad which gives the coordinate transformation between Newtonian gauge 
and the local inertial frame at rest with the baryons. To first order, only e°o = (1 — \l/)/a 
matters. Here ot is the Thomson cross section, and cl^a is the radiation constant. Note 
in the above only the monopole radiation temperature Tr enters at first order. The above 
collision term is nothing else but the Kompaneets collision term for Compton scattering in a 
local inertial frame at rest with the baryon fluid. 



Combining the expression for the energy-momentum tensor (41), (42); the number con- 



servation for each nucleus (36); the evolution of the electron density (37), (38); and the 



component of (47) together with (48) we obtain an equation for the matter temperature, T, 



M- 



To zeroth order we obtain 



M ' zril M 



(o) 



T 



(o) 



a 



M 



' n 



_ A (o) 
(o) A c 



(50) 



The term including Q corresponds to the evolution of the molecular weight (the number of 
particles in the gas changes due to the electron recombination). It has a timescale > 10 35 
larger than the Kompaneets timescale (shown in Fig. [TJ of the term on the right hand side, 
and has therefore a very small effect. The first order expression for the matter temperature 
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is given by 



8T M + 2H5T M + ^T M 6 b + -^Q^5T M 



ru 



+ 



2 a 



3 n {0) V 2 



\T M Q (0) - A 



2t£ } $ + 



n e o e + n n Ob 



n 



(o) 



" " f 5A c - 3 -T M 5Q 



S (°) 



(51) 



To gain some intuition, notice that the 9 b and 7i terms can be combined to obtain the 
local velocity divergence, and that \I/ converts coordinate time to local inertial frame time. 
Neglecting the evolution of the molecular weight, the above equation for 5Tm reduces to 
equation (B.12) in [27]. 



3.3 Discussion of numerical results 



The evolution of S e and 5T M is given by eq.s (36), (38), (40) and (51), where we also include 



the contribution from He recombination as explained at the end of sec. 3.1 The linear 
approximation is valid up to the reionization epoch, when the electron density perturbation 
becomes non-linear. 

The cosmological parameters we use are 

(tt b , tt A , h, T cmb ,F p , n s ,r) = (0.041, 0.76, 0.73, 2.726, 0.24, 0.958, 0.092) . 

We remind that the peak of the visibility function is at 77 = 288.42 Mpc and the present 
conformal time is t)q = 14554.28 Mpc. We find it helpful to give the following correspondence 

r] = (70, 288, 500, 1000, 5000) Mpc z = (5936, 1090, 507, 174, 9) . 

In the Introduction we argued that the amplitude of 5 e at large scales should be enhanced 
due to the small timescale of recombination. This is valid for gauges, as Newtonian gauge for 
example, where the perturbations do not vanish for k — > 0, and are completely described as 
time delayed homogeneous solutions. Therefore, for k — > in the Newtonian gauge we expect 

5 e /5b = d\nn e /d\nnb = 1 — ^(lnx e )/(37^) , (52) 



and analogously o~t m /o~t r = d In Tujd In Tr. In Fig. [3] we show a plot of the amplitudes of 
several perturbations including 5 e and St m in the Newtonian gauge for k = 10 _4 Mpc _1 . The 
electron density and matter temperature perturbations indeed follow precisely an evolution 



equivalent to a time delayed homogeneous universe solution as in eq. (52) (in fact in the plot 
the difference is in practice indistinguishable). The electron density perturbation departs 
from 5b during Helll^ Hell recombination at rj ~ 70 Mpc, and then during Hydrogen recom- 
bination. The Hell— > Hel recombination is delayed [29], and its contribution to 5 e overlaps 
with that of Hydrogen recombination. 
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As expected, the electron density perturbation is enhanced with respect to the other 
perturbations at recombination. The plotted first order variables are (except for St m ) the 
matter and radiation fluctuations which source the second order CMB anisotropies. Later we 
will see that this enhancement persists also for modes well inside the horizon. This suggests 
that 5 e will play an important role in sourcing the bispectrum of the CMB. The calculation 
of these effects will be the subject of a companion paper [28J. 

Still in Fig. [3j before z ~ 500, we notice that St m — St r as Compton scattering keeps 
baryons and photons in thermal contact. St m /$t r becomes equal to 2 after z ~ 500 when the 
matter temperature decouples. At this points matter and radiation redshift differently even 
in an unperturbed universe, which explains the factor of 2. 
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Figure 3: The amplitude of the perturbations to the electron density in the Newtonian gauge 
compared to other first order perturbations for a mode well outside the horizon. The peak 
of the CMB visibility function is at r\ — 288 Mpc. The normalization of the perturbations 
is the one in CMBFAST - superhorizon C, = 1. The perturbations in the Newtonian gauge 
for superhorizon modes are equivalent to time shifted zeroth order evolution. This results in 
two regions of enhancement of S e corresponding to Helll— *>HeII, and the combined Hell^Hel 
and HII— »HI recombination. Similarly, 5t m coincides with St r while Compton scattering is 
effective in keeping the photons and electrons in thermal contant. As Tm decouples, 5t m 
converges to 25t r - 



Even more interestingly, we find that the enhancement to 6 e persists even for modes inside 
the horizon where the effect is clearly not gauge dependent. In Fig. [4] we show 5 e and St m 
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for a mode corresponding to the second acoustic peak. The Kompaneets timescale n t T M /A c 
is still smaller than one comoving Mpc for 77 < 400 Mpc (see Fig. [T|, therefore, for the modes 
we consider, even after recombination, the full (zeroth+first order) Tm and Tr are kept in 
contact by Compton scattering, as can be seen in Fig. |4j At later times, the two temperature 
fluctuations deviate and 5t m becomes greatly enhanced with respect to 5t r since overdense 
regions become hotter due to adiabatic contraction. 

For z < 200 (77 > 1000 Mpc) one can neglect d ri log(x e )/(3H) < 0.1 in (52), and from the 
time shift perspective one would expect 5 e /5b ~ 1. However, from Fig. [4] we clearly see that 
5 e is suppressed compared to 5b during the dark ages. As argued in [26], during the dark 
ages, overdense regions have enhanced T M , and this is associated, as we will soon explain, to 
a more efficient recombination and a suppressed 5 e . This effect induces a decrease of about 
10% in 5t m during the dark ages, which in turn influences the 21 cm power spectrum at a 
2% level at z ~ 50 [26J. For a detailed discussion of the effect of the matter temperature 
perturbations on the 21cm angular-power spectrum we refer the reader to [27]. 
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Figure 4: The amplitude of the perturbations to the electron density in the Newtonian gauge 
compared to other first order perturbations for a mode corresponding approximately to the 
second acoustic peak. The enhancement of 5 e is still well approximated by eq. (52). This 



mode is well inside the horizon, and therefore, this enhancement is not a gauge artefact, and 
will leave an imprint on the CMB bispectrum. 



From the numerical results (see Figures [4] and [5|, we see that eq. (52) is still a good 
approximation to the ratio of 5 e and 5b even for modes well inside the horizon. This does 
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not mean anymore that 5 e is just a time delay of the homogeneous solution. Instead, it 
tells us that n e is determined by the local value of (or equivalently of the temperature) 
with the same relationship as in an unperturbed universe: 5n e ~ 5nbn e /fib- This can be 
explained by the fact that the timescale for recombination is still much faster than the mode 
we are considering. In Fig. [5] we plot 5 e and 5b (which has a size comparable to the other 
perturbations that source the second order CMB anisotropies) as a function of k at the 
peak of the visibility function. As one can expect, (52) well approximates 5 e down to scales 
comparable to the recombination timescale (k ~ 0.1 Mpc -1 ). The enhancement of 5 e is 
prominent for modes with scale larger than approximately 1CT 1 Mpc" 1 (which is comparable 
to the photon diffusion scale kp ~ 0.15 Mpc" 1 ). The suppression at high-/c with respect to 
what implied by (52) is due to the fact that for these modes the timescale of the oscillation 
becomes faster than the rate of recombination (see Fig. |2|, inducing an averaging out of the 
enhancement 

Notice that in Fig. [4j we can see that at late times, even for a relatively low-k mode, 
\5 e \ is a bit suppressed with respect to what is given by eq. (52). Let us briefly explain 
why this happens. For rj > 400 Mpc, the timescale of the mode with k = 0.04 Mpc -1 
becomes comparable to the recombination timescale (see again Fig. [2]). At this time, the 
contribution from Lya escape, two-photon absorption, and photoionization is negligible, and 
we can approximate the homogeneous integrated Hydrogen recombination collision term as 
Q ~ — a^nnn^e. This clearly increases in magnitude for higher n e or rimi- This dependence 
of Q oc — 7Zhii implies that overdense regions will recombine faster, explaining the suppression 
of \5 e \ relative to (52) at later times in Fig. |4j 

Let us comment on the validity of our calculation at high fc's. In deriving the perturbed 
Peebles equation, we pointed out that for k > IMpc -1 , the Lya line cannot be treated as 
quasi-static, and the perturbations to the escape probability will be modified. However, those 
perturbations have a comparably small effect on 5 e , since recombination is dominated by two- 
photon decay. The most important effect which may influence the 3-level atom calculation is 
the effect of nonequilibrium excited levels. That effect in the homogeneous universe leads to 
a correction of about 10% to n e for i] > 360 Mpc, when photoionization becomes negligible 
(Fig. [2]) [29J. At this point recombination to excited states becomes enhanced, and the bound- 
bound transitions between n > 2 levels is no more in thermal equilibrium (for a discussion of 
the effect see [15]). This tells us that the timescale over which the excited states equilibrate 
becomes comparable to the Hubble time at these redshifts (see again Fig. [2]). Since this effect 
occurs quite late in the recombination history, it gives only a 10% effect on n e . In the presence 
of perturbations with k lower than ~ 10~ 3 Mpc" 1 , the effect will be of the same order. For 
modes with higher k, the thermalization timescale becomes more important even earlier, and 
so potentially the effect could be larger. However, since the photoionization rate grows quite 
steeply as we move to higher redshifts (see again Fig. [2|, this effect becomes important only 
for k ~ 1 Mpc" 1 . This is the upper limit on the /c's for which our calculation is valid. 



9 This suppression at large k will lead to an important simplification when we derive the CMB bispectrum 
from recombination 281. 
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Figure 5: A plot of some first order quantities in the Newtonian gauge at the peak of the 
visibility function showing the enhancement of 5 e . The first order 5^ shows the typical mag- 
nitude of the variables sourcing the CMB anisotropics at second order, and one can see that 



8 e dominates. We also plot the expected electron density perturbation coming from eq. (52) 
(dotted line) which is a good approximation to the true S e down to the recombination scale 
(k ~ 0.1 Mpc -1 ). The enhancement for superhorizon modes is gauge dependent, but per- 
sists for modes inside the horizon, where the effect is observable, approximately down to the 
photon diffusion scale (hp ~ 0.15 Mpc -1 ). 

3.3.1 Perturbations in the synchronous gauge 

In [28] we will be working in synchronous gauge to calculate the CMB bispectrum generated 



from recombination. In this gauge, the evolution of S e and ST M is given again by (36), (38) 



(40) and (51), but with the substitutions ^ — > 0, and $ — » —h/Q, where h is the trace of the 
perturbations to the spatial part of the metric as defined in |44J. 

In Fig. [6] we show the results for 5 e in the synchronous gauge. As expected, for modes 
well inside the horizon, the density perturbations in the synchronous and Newtonian gauges 
are nearly identical (cf. Fig. [5]). Perturbations vanish for superhorizon modes since they 
correspond to time-shifted homogeneous solutions and the synchronous gauge time coordinate 
equals the proper time of comoving observers which makes the time-delay vanish outside of 
the horizon. As we will show in [28], this property of the synchronous gauge will be very 
useful in estimating the CMB bispectrum generated from recombination. 
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Figure 6: The same as in Fig. [5] repeated in the synchronous gauge. The only important 
difference between the results in the Newtonian and the synchronous gauge is for superhorizon 
modes for which the amplitude in the synchronous gauge goes to zero. 

4 Perturbations to Second Order 

We now begin the second part of this paper, where we will give the full set of second order 
Boltzmann and Einstein equations that are necessary to compute the non-Gaussian signal 
introduced in the CMB by the standard cosmological evolution. In the first part, we computed 
the first order perturbation for 5 e because, as we will see, it will be a source for the second 
order perturbations. In this derivation, we will follow quite closely BMR1 (i.e. Ref. [T§]), 
who first gave the full set of second order equations, though in our derivation there will be a 
few important differences. 

4.1 The Poisson gauge 

The metric in Poisson gauge [H] is given by 

ds 2 = a 2 (r]) [-e 2 % 2 + 2u i dx i dr] + (e~ 2 % + x^dxW] , (53) 

with Xa = 0; since the trace can always be absorbed in The gauge fixing is done by 
choosing 

= , Xij,j = , (54) 
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for each order in perturbation theory [49J. Following [Tj5] we restrict the Poisson gauge by 
choosing the vector and tensor modes given by uJi and Xij to be second order quantities 
(we neglect here primordial vector and tensor modes). For convenience in the expressions 
that follow, we will often use e 2 * and e 2 *. In this gauge, there is one scalar degree of 
freedom eliminated from and one scalar and one transverse vector degrees of freedom 
from gij. Thus, u>i is a transverse vector, while \ij is a traceless-transverse tensor. The scalar 
potentials \I/ and $ contain the first and second order potentials as follows: \I> = + \&( 2 )/2, 
$ = $(i) _|_ $( 2 )/2, where <3> (2) and \l/ (2) are meant to be of order (<3> (1) ) 2 . We notice also the 
nice property that \f—g = e* _3 *a 4 to second order. Notice that, with respect to BMRl's 
conventions, we substituted $ <-> so that they match the usual notation in the Newtonian 
conformal gauge. We use the time-independent piece of the spatial part of the homogeneous 
FRW metric, <5 lJ and 5^, to raise and lower spatial indices. Thus, we do not keep track of the 
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placement of spatial indices 

In order to derive the Boltzmann equation it is useful to find the change of coordinate 
from Poisson gauge to a local inertial frame. This can be found by solving for the tetrad e M a , 
which satisfies 



(55) 



Solving to second order we get: 



UJi 



e\ = 



6 ^ij 2 



(56) 



This tetrad is not unique since the local observer coordinates can be Lorentz transformed. 
The momenta in the two frames are related by the tetrad in the following way: 



v'P ' 



(57) 



where P M are the momenta in the Poisson gauge, and p^' are the ones in the local inertial 
frame. The tetrad above is chosen such that P l = for particles at rest with respect to the 
inertial observer (i.e. when p — 0), which allows for an easy interpretation of the observers' 
trajectories (see below). Using (56) and writing p^p 1 = p 2 and p l = pri 1 (such that n l rii = 1), 



we obtain the following relation between Poisson gauge momenta and local inertial frame 



10 Since in this section we are summarizing second order results in Poisson gauge, we find it useful to mention 
that BMR give the correct second order Christoffel symbols for this gauge in eq. (A. 2) in their [21] (BMR1 
have instead a sign typo in their eq. (2.8) for the Christoffel symbols). 
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momenta 
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p o 



—e~ x 
a{rj) 



P_ 
E 



LOiU 



P eV I 6. 



Pi = ae [pi + EuJi + -XijP 3 



(5f 



The observer's trajectory is given by dx 1 / dr]\p = Q = P l /P°\^=o = 0. Thus, the local inertial 
frame that we have defined in this way is nothing but the Local Inertial Frame Instantaneously 
at Rest with respect to Comoving Observer (LIFIRCO). 



4.2 Streaming terms 



We are now ready to begin to write down the Boltzmann equation. From ( 110 ) of App. A we 



can write the Boltzmann equation for the diffeomorphism invariant one-particle distribution 
function, /, as 
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yQ dJ =p ^dl + dP^dl 

drj dxi 1 dX dP % 



C[f] 



(59) 



where C[f] is the diffeomorphism invariant collision functional. As we explain in App. [A] 
and as we have already used in the former section, it is useful to introduce 'nice' coordinates 
where the momenta are the ones of the local inertial frame, and we define 



F{x*,it) = f{x»,P j (p\x»)) . 



(60) 



In this way the collision term is the same as in Minkowski space, and all the metric pertur- 
bations are confined to the free-streaming term. In these nice coordinates, the Boltzmann 
equation becomes 



dF dF dx i dF dp 

1 : 1 r — 

drj dx l drj dp 1 drj 



1J L = ^C[F] 



or equivalently 



13 



dF dFdx i OF dp dF dn l 

+ 7^"^ + ~ 7- + 



(61) 



(62) 



dr] dx l drj dp drj dn % drj P° 

11 Note that the expression eq. (3.6) for P l in [T5] is not a linear function of p M . Thus, their ?/"s that 
are present in the streaming terms are not momenta measured in a local inertial frame. This leads to an 
inconsistency since the p M appearing in their collision term should be the momenta measured in a local 
inertial frame. Thus, the streaming terms derived in the next section are different from the ones in BMR1. 
However, it turns out that this has no effect on the streaming term of the final equations for the evolution 
of the baryon and CDM number density on their continuity equations, as well as on the photon brightness 
equation at second order, because of what seems to us to be just accidental cancellations. 

12 In order to help comparison with the existing literature, we notice that this equation is equivalent to 
eq. (4.1) in |21j . but differs from eq. (3.1) in [TS]. 

13 For simplicity, we denote both F(x\p';r]) and F(x > ,p,n°;ij) with the same symbol F. It is clear which 
one is used in all equations below. 
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To derive the Boltzmann equation to second order, we need to calculate dp/dn, dp 1 /drj 
and dx l /dr] to second order, and dri /drj to first order, since OF /dri is already a first-order 
quantity. 

Writing the coordinate velocity dx l /drj to second order is simple: 



dx % 
dr\ 



pi 



— n J e 



<h , ( 1 - ^ k n k J - 



(63) 



The ijj piece in (63 ) converts E as observed by LIFIRCO observers to the energy of the particle 
as measured in the normal inertial frame, defined as the frame which is free falling, which 
moves with a velocity (— uj 1 ) with respect to the LIFIRCO observers [IS]. The rest of the 
terms appearing above convert from inertial frame velocity to coordinate velocity. 
To obtain dp 1 /drj and dn % /drj we write the geodesic equation 



dPv 
drj 



pu pp 
vp pO 



which, by using (58), and after some simple manipulations, give: 

' E 



dri 1 
drj 
1 dp 
p dt] 



- (5 lj - 
---H + <5> 



rin j ) 



P 



E 



~-n l n 3 Xij 



E 



m 



-n uoj 



-UuJin 1 



(64) 

(65) 
(66) 



p 2 p Ep 

As one could have expected, dri 1 /drj is affected only by the transverse gradients of the scalar 
potentials. The second term in (65), which for nonrelativistic particles is 0(v 2 ) suppressed 
compared to the ^ j term, is responsible for the fact that relativistic particles are deflected 
twice as much compared to nonrelativistic particles. 
We can also write dp k /di] to second order 

1 dp k 
p drj 

E. 1 t. P 



-n k H + n k § - e 



p E 



^n h ri$ ti 
E 



m 
pE 



Hu) k + 



+n [ujih - 0Jkj, 



2E 



n l n> (xik,j ~ Xij,k) 



p 2 

where a dot denotes 8/ drj, and H = a/ap^j 

14 For the sake of comparison, we can rewrite the above equation in a form similar to equation (4.70) of 



(67) 



drj 



aE 



1 + \x- ! v 



& + -^P x (V x lj) 



P\ 2 



E 



d n {aE) uj 



(68) 



where we defined [x ■ p\i = XijP 3 ■ The above equation is correct to second order in the metric perturbations. 
The reader is referred to equation (4.70) of [H] and its discussion for the interpretation of the terms appearing 
above. Since here LOi and Xij are second order, and ^ and $ contain the second order scalar perturbations, 
the only additional component to (4.70) of [48] is the extra factor of e~* in (68 1 which converts coordinate 
time to proper time. 
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4.3 Evolution of the number density and momentum conservation 
at second order 



The particle 4-current measured in the LIFIRCO, which we call iV M , can be written using 



eq. (112) of App. A 



N» = g dcg 



d 3 p 



The average fluid velocity v 1 measured in the LIFIRCO is given by 

^wo- 



rn 



(70) 



It is good to know that the transformation from the LIFIRCO to the local inertial frame 
momentarily at rest with the fluid is given by a boost with velocity v i q| This implies the 
standard relation N l = N°v % = n v 1 / yl — v 2 . To get an equation for N°, we can follow BMR1 



and multiply both sides of (61) by P°/E and integrate over d 3 p. Another approach which 



simplifies the calculation is to use ( 109[ ) and (113) of App. |A]to write the conservation of the 
particle 4-current as A/" M M = f it C[f] which, by expressing A/" M in Poisson gauge in terms of 
by the tetrad relation (156]), to second order gives: 



9d. 



N° + 3N°(H - $) + e*+*N} + e*+*JV < (* )< - 2$,,) 
d 3 p 



E 



-C[F] . 



(73) 



In this equation all the variables are meant to be the full variables, containing their zeroth, 
first and second order pieces 16 Notice that for the perturbation 5b to the nuclei number as 



defined in (43) the right hand side above vanishes. 



Next we derive the momentum conservation equation for baryons (meant as usual as the 
collection of ions, electrons, and atoms). In doing that, one has to be careful in treating the 

15 Under a boost with velocity v we can go from the LIFIRCO to the instantaneous rest frame of the fluid. 
Let us check that the average velocity v of the fluid coincides with v. We can use the invariance of d 3 p/E 



( 102 1 and / to rewrite (70 1 as 



fdcg 



d 3 p' 

—f(x,P j ($,i?;r));r))p l 



(71) 



Using the transformation of p M under a boost with velocity v such that we go to the rest frame of the fluid, 
we can express p in terms of p 1 and E' , with which we denote the momentum and energy of the particles 
measured in the fluid rest frame. Assuming that the distribution function has vanishing dipole in the rest 



frame of the fluid (by definition of the fluid rest frame), from (69 1 evaluated in the fluid rest frame and (71 1 
we obtain 

F=v. (72) 



Thus, v coincides with the average velocity of the fluid v. 

16 This equation is also given by BMR1, but they set the right hand side to zero even for n e , which is an 
approximation that during recombination is not satisfied. 
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momentum transfers in radiative transitions that do not affect relevantly the kinetic energy 
of the fluid, and in the Coulomb interactions, which force electrons and ions to have the same 
velocity. We already discussed these effects when we derived the perturbations to the kinetic 



matter temperature. Thus, we derive momentum conservation from the ith component of (47) 



which involves only the photon second order Compton collision term (given in our eq. (139) in 
App. [C]). Unless stated otherwise, we denote the order in perturbation theory of all variables 
by a superscript in parenthesis, which we drop whenever this does not cause any confusion. 
Also, for simplicity, we drop the subscript & from the baryon velocity Thus, the first and 
second order expressions for the baryon momentum are given by: 

+ 9?' = ^ („(«. _„<*) , (74) 

which coincides with BMR1, and p] 

1 / ( q„,(2)i \ 

- ( + HvV* + 2u/ + W.J + - ¥%M + v^v W i + 

2 \ drj ' J ' 3 



A n* ' r / \ 

+(*« + d>(D)*« = If ^ [(a« + * + S e - 8 b ) („«* - «««)+ 

,(2)< ,,(2)< 



(0) 



2 2/4 



(75) 



In analogy to BMR1, we have found it useful to define the photon velocity as 

1 

P-y+Pj 



< = I / (PpFyP* ■ (76) 



f = — n e cttcl is the differential optical depth, Aq = 5~ p and UPj' is the photon quadrupole 
defined as 



(77) 



We also neglect terms very small suppressed by Tjv//m p , where m p is the proton mass. The 
term proportional to Aq + 5 e — 5& on the right hand side is the perturbation to fp^ / , and 
as usual \P converts coordinate to proper time. The velocity measured by comoving LIFIRCO 
observers appearing on the left hand side is converted to the velocity in the normal frame, 
v+u. The combination v+v l v^ converts the Eulerian time derivative to Lagrangian derivative. 
The second order piece of e + ^ corresponds to the contribution to the acceleration from 
the gradient of the Newtonian potential with the appropriate conversions between inertial 
and comoving coordinates. The $ term accounts for the usual redshifting of momentum in a 
time- varying potential. 

17 Note that the left hand side of this expression is the same as eq. (6.44) in [19] (with the substitutions 
\& — > $ — > W). However, the right hand side has acquired a dependence on 6 e . 
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4.4 Photon brightness equation 

From the Boltzmann equation for the photons, we get the second order brightness equation 



^[A»+«»] + £ [ A (D + «W] 



(2) 



4AW($W-$ ( ! 1 V) 



(7f 
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where we have used the Compton collision term given in eq. (139) and we have explicitly 



reminded that in the second term in the first line one has to take just the second order 
contribution. Here we have defined 



x, 77, n) 



j dpp 3 F^ 



fdpp 3 F(°) ' 

where we have expanded the distribution function as 

F(x,p\ri) = F^{x,p\ri) + F^\S,p\-q) + ^F^(x,p\r]) 



(79) 



JO) 



where F^ is the first order perturbation and F^ is the second order one. The p l 's are the 
the moments as measured in the LIFIRCO. The expansion in sperical harmonics of a quantity 
A(n) is given by: 



A 



1 n 1 



-1 



4.5 Energy-Momentum tensor 

We now give an expression for the second order energy momentum tensor and for the Einstein 
equations in the remainder of this section. In what follows, primed quantities are meant to be 
evaluated in the particular inertial frame that is at rest with the fluid. The 00 component of 
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the energy- momentum tensor given in (44) can be written for the matter components using 



(58) in the form 



o - -^de g / d pEb - -g dcg j —E b - -g dcg j — — — % b - - - _ ^ , 

(82) 

where p and p are the energy density and the pressure of the fluid. For the photons we have 

T° o = -2 j d 3 pEF 1 = -pf (l + A« + i Ag>) . (83) 

The other components of T^ v are given by 

T = -g Acg e^ j d 3 P Fp\ (84) 

lp 2 



d 3 p 



p l + u l (E + 



3E 



E 



■Ep l pj . 



5) 
6) 



Let us elaborate on these expressions. From the last equation above we can extract the 
expression for the pressure p: 



v 3 « 



iij,=0 frame 



1 fd 3 p' ajp, 



57) 



confirming what we used in (|82|). For massless particles, this automatically implies the usual 

(88) 



p = 3p. Using (87), we can write to second order 

J d?pF P i + t o i (p + p) = _ e -2(*+*)r o + j( p+ p) 



rpO _ „ 



9de g e 



For non-relativistic particles, the above integrals can be evaluated in a variety of ways. One 
can work with the distribution function in the rest frame of the fluid, following the steps 



we used to derive T° in (82). Or one can expand the relativistic Boltzmann distribution 
function in the comoving frame. Yet another way is to write the energy-momentum tensor in 
the standard covariant form 

rt = (P + P)tW, + ^P, (89) 



and then use (58) to write down the 4- velocity to second order in the metric perturbations: 

U° = a- 1 e~*(l + ViU* + v 2 /2) , U' 1 = a~W . (90) 

This last method is easiest to work with, and by using it for the baryons and the dark matter 
we obtain: 



r = -e*+*(p + P K, 

TV = 5}p + (p + p)vV . 



(91) 
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The above T l is correct also for the photons upon replacement of v with tL. The result is 



different from the one in [T5]. From (86), we recover the result of [IH] for T„ 
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4.6 The second order Einstein equation 

We now give the Einstein tensor, G^, to second order in the Poisson gauge: 
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(93) 
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(96) 



with the projection operator defined as P^ = [SikSji— 6^6^/3], and k 2 = 8tcGn. The equations 
above are equivalent to the ones in App. A of J2U] with <3>here = ^thcrc and ^hcre = $ there- 
The transverse part of the last equation above gives the gravitational waves generated by the 
scalar mode [50] . 

The remaining longitudinal part of G l j = kT^ is given by 

Q , 

where Q is defined through the following chain of definitions: 

V 2 Q = 3N-P , 

with 

72 ■ 



P = P\, and V z N = P* hl3 , 



where finally: 



Note that our expression for P*. is equivalent to the one in 
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(97) 
(98) 
(99) 

(100) 



5 Summary 



The purpose of this paper was to provide a detailed analysis of the perturbations to the 
recombination history of the universe and give the full corrected second order Boltzmann and 
Einstein equations for the photons, baryons and CDM for a flat FRW background. Those are 
necessary for the correct treatment of CMB secondaries coming from recombination. 

We analyzed the different timescales in recombination and argued that a perturbed version 
of the recombination equation for the Peebles 3-level atom is adequate for obtaining the 
evolution of the perturbations to the electron density. We derived rigorously the perturbations 
to the escape probability by angle averaging the photon distribution function, and showed 
that it is sufficient to treat the escaping of Lya photons as only due to the local velocity 
divergence. 

We showed that for modes longer than the timescale of recombination (corresponding to 
k < 0.1 Mpc -1 ), it is a good approximation to consider ne° +1 ^ as determined by the local value 
of ^ (or equivalently of the temperature), with the same relationship as in an unperturbed 
universe: 5n e ~ 8n\,h e jh^. This is due to the fact that the timescale of recombination is much 
faster than the considered modes. For the same reason, this also implies that the amplitude 
of the electron density perturbations is enhanced by a factor of roughly five compared to the 
baryon density perturbations. This enhancement is physical since it persists for modes well 
inside the horizon, down to scales comparable with the photon diffusion scale (corresponding 
to k < 0.2 Mpc" 1 ). The enhancement is most prominent around the peak of the photon 
visibility function and therefore has the potential to have an impact on the CMB bispectrum, 
as we will investigate in a companion paper [28] . 

Solving for the density of free electrons in the universe, as well as writing the full set of 
second order Boltmann and Einstein equations requires to be careful with the definition of 
coordinates and collision terms. For this reason, we reviewed the formulation of kinetic theory 
in curved spacetime and applied it to derive the Boltzmann and Einstein equations at second 
order in the perturbed FRW universe. We compared our results with former pre-existing 
literature. We find that the former derivations are not completely correct, and we give the 
necessary corrections. 

Solving the full set of second order equations is a very hard task, which however is mo- 
tivated by the current experimental and theoretical status. In a companion paper [28] . we 
solve approximately these second order equations restricting ourself to modes well inside the 
horizon and to second order sources that are proportional to 5 e , and we compute the induced 
bispectrum. There we find that the non-Guassian signal is equivalent to about /nl ~ — 5, 
which is potentially detectable by future experiments such as Planck. 

Note Added 

While this paper and its companion [2B] were slowly written down, preprint [5T] appeared 
very recently. Our results about the full set of second order equations in sec. [4] partially 
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overlap with those of [5*T] . 
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A Covariant Form of the Boltzmann Equation 



A.l Preliminaries 

In this Appendix we intorduce some useful quantities in covariant kinetic theory following 
closely the work of Ehlers [46] , to which we will often refer for details. We work in generic 
coordinates, and we define the affine parameter of each particle trajectory so that 



~d\ 



(101) 



where P M is the physical momentum. In practice A of one particle is defined so that dX 
coincides with the time interval of a local inertial observer divided by the energy observed in 
the same frame: d\ = Gfa inertia i/P inertial . Let us first introduce the diff. invariant measure of 
the one-particle phase space. The measure in momentum phase space is given by: 



vr = g Acg 2^d A P l J dP u 5(P^ - m 2 )0(P u 
2d 3 P 



9de 



9 d 3 P l 



2cr P C 

9de S ^= J dPoSiPxp^ - m 2 )6(P ) = g deg - 



\Po\ 
d 3 P t 

P°\V=9 



(102) 



where c^eg are the internal degrees of freedom (degeneracy) of the particle species for which 7r 
is calculated. The relationship between the momenta in generic coordinates and the momenta 
in the local inertial frame is given by: 



(103) 



Here, p° = —po = E is the instantaneous energy of the particle with respect to the local 
observer. The tetrad satisfies the set of equations e^e^,^ = r}a>pi. For example, in the par- 
ticular case in which the generic coordinates are the ones of Newtonian gauge, they are given at 



second order in eq. (56). From (103), we can rewrite dP l = [e l Q ^-\-e % 3 \dp^ using EdE = pjdpi . 



Vwa' (e 



p v , it is 



Plugging this expression for dP l in (102) and using P 
straightforward to show that 

d 3 p 

We find that, once expressed in terms of momenta in the inertial frame, the momentum-space 
measure reduces to the standard Lorentz invariant one. The one-particle phase space needs 



(104) 
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also the measure for the spacetime volume. This is given by the standard 4-form w = ^/—gd^x, 
where g = detg^. The one-particle phase space is 7-dimensional, and its measure is given by 
Q = wk. Notice that both the spacetime volume measure and the momentum-space measure 
are separately diff. invariant. 

A. 2 Invariant distribution function and the Boltzmann equation 

The Liouville operator L is a vector acting on the 7-dimensional space (x^, Pi). It is given by 

d 8 dP l 8 

L clearly describes the free streaming of a particle in the 7-dimensional phase space. If we 
contract L with Q we obtain a 6-form: 

uj = L ■ Q . 

Up to a constant factor, which will be fixed shortly, this is the unique 6-form which assigns 
a nonzero volume to any hypersurface not tangent to the vector flow L. By construction, 
to is invariant with respect to the phase flow, i.e. J s uj = L, o>, where E and E' are two 
cross-sections of a tube, defined by the Liouville flow. Doing the contraction in the expression 
for uj one obtains 

uj = p^ ^ , (106) 



where a a is the covariant surface element 18 For our physical interests, we restrict uj to be 
defined on surfaces £ of the kind £ = G x K, where G is a spacelike hypersurface, and K is the 
mass shell hyperboloid in momentum space. Going to a local inertial frame instantaneously 
at rest with a particle with 4-momentum P M , we see that uj reduces to d 3 xd 3 p. Thus uj is the 
standard 6-dimensional phase-space volume measure, a fact that fixes the remaining constant 
that made the definition of uj not unique. 

If we call iV[S] the average number of 7-dimensional trajectories that cross a certain surface 
E, we can define the phase-space one-particle distribution f(x\ Pf, rf) through the relationship 

N[E] = [ fuj . (107) 

By varying E and making it arbitrarily small, and by using the uniqueness of uj, one can show 
that / is a uniquely defined scalar function jl6] . 

The average number of interactions in a phase space region D, N[dD], is given by 



N[8D] = [ fuj= [ d(fuj) = [ L[m ■ (108) 

JdD JD JD 



18 Given a set of coordinates (u 1 ,u 2 ,u 3 ) parametrizing a 3-surface, the 3-surface measure can be written as 
o-q = {I /2>^^ge a( 3 1 s{dx 13 /du 1 )(dx' y /du 2 )(dx 5 /du 3 )d 3 u with e a p jS = £[ Q /3 7 5], and e i 2 3 = 1 [52]. 
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In the last passage we have used that d(fu) = L[f]fl (see [IS])- Thus, the phase-space density 
of collisions is given by L[f] = df /d\. With this interpretation of L[f], we can finally write 
the Boltzmann equation in covariant form: 

L[f] = C[f] , (109) 

where C[f] is the collision term giving the phase-space density of collisions, which we will 
discuss in the next subsection. Before moving on, let us notice that in the above equation 
one can use L[f] = df/dX = P°df /drj to obtain 

J-|-ow, (no) 

which is the form of the Boltzmann equation which we use in most of the paper. 



Notice that we can rewrite ( 108 ) as: 



A' 



N[dD]= a, / ttP"/ , (111) 
Jg J k 

which allows us to read off the expression for the particle current as 

JV> = / irP^f . (112) 
We notice that its divergence is given by 

M%= I _L[f]n. (113) 

A. 3 The collision term 



A 



Having introduced the invariant distribution function, /, and the Boltzmann equation (110), 
we can procede with the discussion of the collision term, C[f], which is a scalar by construc- 
tion. Let us consider interactions of the kind iKi,jK 2 ■■■—>■ rK 3 . . . , where ■ ■ • , r, s, . . . } 
label the species of the particles, and K a labels the momentum regions for each particle. The 
collision term for /, is then given by 



C[fi] = [ ny f 7r r ---[f l (x k ,P 1 ; V )f j (x k ,P 2 ; V )---(l±f r (x k ,P 3 ; V ))---] 
Jk 2 Jk 3 

xW(g^iP u jP 2 ,--- ;riV") , (114) 

where P n is the 4-momentum of the n th particle participating in the collision. Notice that in 
the equation above there is obviously no integration over 7Tj. Since the phase space measures 
are diff. invariant, so it is W(g^ u ; {P})- We can write it in a local inertial frame W(g^ u ; {P}) = 



W(t] i1u ; {p}), where the transformation between p and P is given by the tetrad in eq. (103). 
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In a local inertial frame W^^p) is expressible using the usual scattering operator S: 

W(r) IMU ;ip 1 ,jp 2 ,--- ;rp 3 ,---) = (27r) 3m - 4 5 (4) (p 1 + p 2 + p 3 ) 

x \ (rp 3 ,--- \M\ip 1 ,jp 2} ---}\ 2 , 

(rps, • • • \S - 1 \ipi, jp 2 , ■ ■ •) = 5 (4) (pi+p 2 H P3 ) 

x (rpa,--- \M\ip u jp 2 ,---) , (115) 

where m is the total number of particles in the interaction, and Pi is the momentum of the 
i th particle. 

A. 4 Nice coordinates 

Notice that the collision term is a scalar local quantity in spacetime, and therefore, if we write 
it in the frame where the momenta are the inertial ones, then it does not contain any metric 
perturbation. This is a nice simplification that we wish to exploit. In order to do this, we 
can define a new distribution function F(x^,p l ) as: 

F(x M , p 1 ) = /(^,PV,«) • (H6) 
The new distribution function is defined everywhere in spacetime, and in each point it is a 



function of the local inertial frame momenta defined by eq. (103). We call this 'mixed' set of 



phase space coordinates 'nice coordinates' 19 In fact, in these coordinates the collision term 
C[f] becomes just the standard Minkowski one, with no metric fluctuations. In particular, 
the momenta measure just becomes 

7T; = £deg,i-^ , (117) 
and the metric which contracts the momenta becomes rj^. Coming now to the Liouville term 

rf/(^,py,^) 

Tx ' (118) 

it is easy to write it in terms of F{x^ 1 p l ) by applying the chain rule. Since the Liouville term 
is a total derivative with respect to A, this simply amounts in doing 

rf/(x^,P i ) dF{x^,p i ) 



dX dX 
The Boltzmann equation in nice coordinates becomes: 



(119) 



dF dFdx^ dFdf = l 

dr] dx i dr] dp* drj P° 1 1 ' 1 ' 

which is the version we use in the paper. All the dependence on the metric is now just inside 
dx l /dr), dp 1 /drj and P°. 



19 This is what is usually done, in a less formal way, when studying the Boltzmann equation at first order 
|53j . The generalization to second order can therefore be a source of confusion. Here we are explaining how 
to perform correctly this change of variables directly at general non-linear level. 
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B Derivation of equation (40) 



B.l Einstein coefficients 

We are now going to derive the expression for the escape probability at first order in the 
fluctuations. In order to do this, it is useful to first derive an expression for the Einstein 
coefficients using the notation we introduced in the previous sections. For absorption of a 
photon by an atom in a transition from the i th to the j th level, we have the following spacetime 
density of interactions 



iMlAbs. 



^j^jf-yfiW(i,r,j) 



(121) 



where f\ denotes the one-particle distribution function of the atoms in the i th level, and 
/ 7 is the one-particle distribution function for the photons. Going to the inertial frame 
instantaneously at rest with respect to the the atoms (£f 6 = 0), i.e. to the laboratory frame, 
the result for the spacetime density of absorptions is given by the following term of the 
standard radiative transfer equation 



>M I Abs. 



B{jTl{ 



B- 

I(v,h)(j)ij(v,n)dvdri/4:7r = g~. — — rii I vF^iAv^n)d i 'p 

An 



(122) 



where g 7 = 2 is the degeneracy of the photons, Bij is the Einstein coefficient for absorption, 
v and p are respectively the energy and momentum of the photons measured by an inertial 
observer for whom the bulk velocity of the atoms is zero, <pij{v, h) is the laboratory absorption 
line profile, satisfying J (j)ij{y,n)dvdVt/An = 1, and n« is the proper number density of atoms 



in the i th level. In (122) we used I(u,h) = g 7 z/ 3 F 7 , where / is the light intensity. 



Equating (121) and (122), we obtain an expression for B^ in terms of a collision term: 



^(u,n)B ljni 



KiTTjfiW(i,r,j] 



9i9j 



d 3 pi d 3 pj 
Ei Ej 



FiW(i, r ,j) 



vj,=0 frame J J ~" i ^3 i%=0 frame 

(123) 

where and gj are the degeneracies of the i th and j th levels, respectively. The integral over 
7r 7 can be dropped because this equation is valid for any / 7 . 

The derivation of the Aji and Bji Einstein coefficients proceeds in an analogous way. The 
spacetime number density of downward transitions is given by: 



Af: 



BjiUj / I(v,n)<pji{y,n)dvdVL/A'K + Aj^rij / (pji(u,n)dudil/A7r 



A 4 1 Em. 



9i^ n j \ vF 1 (p ji {y,n)d p+ 



7*3 ji 9 

bjid p = -^-n-iPjA 



A-k ij / YJ1 V 7 B 



But it can be also read off from the Boltzmann equation: 



(124) 



3 j 7T 7 . 



» A* I Em. 



Vf,=0 frame 



(125) 
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Equating these two expressions we end up with an analogous expression for Bji and Ajf 



V)j=0 frame 



= 9i9j 

i?(,=0 frame 



d 3 pi d 3 pj 

Ei En 



A 



.1 1 



Vf,=0 frame 

(126) 



The last expression gives just the standard relation A^ = g^Bji. We can also see that the 
laboratory profiles for stimulated and spontaneous emission are identical, given by <j>ji(v, h) pi 



B.2 Perturbations to the escape probability 



The escape probability is the only non-local interaction during recombination which enters 
in the simplified Peebles treatment, so it needs a special discussion. In the homogeneous 



universe calculation in Section 2J2 we derived the escape probability from the Boltzmann 
equation which we in turn obtained from the detailed balance of the number of photons. In 
this section we will instead start directly with the Boltzmann equation. 



The Boltzmann equation (31) and the expression for the collision term (114) written for 
the photons give 



dX 



^ [/,-(! + / 7 )W(j;i,7) - LfAYii.-jY 



(127) 



Using (123) and (126) we can rewrite the above equation in the inertial frame at rest with 



the baryons to first order as 



dX 



1 



8nu 



[AjiUj - 2v i {B ij n i - BjinAF^ 



:i28) 



f7(,=0 frame 

where A is an affine parameter, and Vj, is the baryon bulk velocity. Note that in a homogeneous 



universe the above equation reduces to (17). 



As in the homogeneous case, we can reparametrize time by using the frequency v(t) of a 
photon as measured by inertial observers at rest with the baryons. Since we are working in 
a perturbed FRW, we need to specify a direction of propagation n for the photons. Thus, all 



quantities in eq. (128) acquire a dependence on h further than the dependence on v. We can 



integrate it to obtain the analogous of eq. (]21|): 

FJu B ,n)e~ T ^ 



F^u, h) 



+ e~ 



VB 



2v 



/3 



dv 



(129) 



which just comes from the integration on each separate direction of eq. (18). In order to 



compute the evolution of the number density of electrons (see for example eq. (122)), we just 
need an expression for the angle averaged escape probability. Since we are working to first 
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This holds in the semi-classical approximation to first order in the intensity 
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order, and all quantities are angle independent at zeroth order, in the Sobolev approximation, 
we find that the angle averaging factorizes as follows 



(F,) n = P(F bb ) n + (l-P) 



p = 



dv(<f>(v,h))nexp[-(T(v,h))n] = (7>n (l-e' 1 /^) 



where 



j(v,h) = 



dv 



(130) 



(131) 



dX u 2 (BijTii — BjiTij) 

Here angular brackets denote angular averaging, and a bar denotes frequency averaging over 

fdv<f>(v)) 



the line profile (i.e. 



This is the analogous of eq. (29) that we found in the 

-Hv 2 



unperturbed case. Notice also that 7 reduces to the homogeneous result when du/dX = 
as in the unperturbed Universe. In the Sobolev approximation the angular dependence in 7 
enters only through the factor dv/dX, and therefore we need to find its angle average to first 
order. Let us rewrite it as 



dv 
dX 



p c 



dr] 



(132) 



Here v is measured in the rest frame of the atoms (we are avoiding the subscript ^ = o for 
clarity), while P° is the momentum measured by a comoving observer sitting at fixed x. If 
we introduce the momentum measured in a LIFIRCO p, the transformation between p and v 



is just a Doppler shift, which to first order is given by p = v(l + v b • n). Combining (58) with 

d(p(l - v b ■ n)) 



(132) we obtain 
dv 



ae*P° 



dv 

dr] drj 



v(l + v b ■ n)- 



drj 



v{\ +v b -n)[ ^(1 -v b -n) 



dv b 
dr] 



nv 



(133) 



Combining this with the geodesic equation in the Newtonian gauge (see e.g. 
first order we obtain 



or (66)), to 



ae*P° 



dv 
dr] 



-Hv 2 



$ .\I> , 

1 - T7 + n % —j + v b ^n l + hi 
ri ri 



dv\ 1 
dr] TC 



(134) 



All first order vector quantities can be written as the gradient of a potential, i.e. A = VA 
Here we are interested in distances £ small compared to the typical gradients of the first order 
quantities, and thus, when angle averaging, we can Taylor expand around the origin: 

A(x) = A(x = 0) + th%A(x = 0) . 

Notice that we are ignoring the time dependence of A. If included, this would give similar 
effects to the ones we will find as due to the space dependence, which as we will see shortly, 
give negligible contribution. When angle averaging n l A;(£n), therefore, we obtain 

(Zn i n?a i d j A{£=0)) a = = ° } . 
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Thus, after angle averaging (134), we find 



d V/ n 



His 2 



$ £V 2 ^ £ € • 1 

x& "^ + 3^r + 3^ + ^ 6 + ^ 



(135) 



where we denoted 9b = djV J b and used dv^/drj = v b + v^^dx 1 jdr\ — {T 6 + raV 6ii (where we 
applied that to zeroth order dx l jdr\ = h\ since Vb is already first order). In (135), because 



of the weighting with the line profile in (130), £ is order Ls- Since £ multiplies first order 
quantities, the above equation holds since to zeroth order we have d^/drj = 1. Comparing 



the above equation with eq. (40) we can see that the perturbation to the escape probability 
has acquired contributions other than the local velocity divergence, Uft/3 = H. Combining 
(135) with the momentum conservation equation for the baryons in Newtonian gauge jH], we 



obtain 



d V/ n 



v 



■c 2 s V 2 S b + \—an e a T {9^ 
& Pb 



(136) 



Thus, the additional terms in (135) correspond to the change of the momentum of the baryon 



fluid due to pressure forces and Thomson scattering in a physical time interval equal to a£. 
The terms multiplied by £ are nonlocal effects due to the photons sampling the Sobolev patch. 

The Sobolev approximation holds only when the perturbation scale k~ l is much larger 
than the Sobolev scale. Comparing numerically the terms in (135) and using £ ~ Ls ~ 5 kpc 



we identify the non-negligible first order terms. We find that within the Sobolev timescale, the 
change in the baryon momentum due to pressure forces and Thomson scattering is negligible. 
Therefore, we can write 



dv\ 
d\/ 



Hv 2 



1 - 




(137) 



Combining (137) and (131) we find that the angle-averaged Sobolev parameter is given by 



(7>i 



(138) 



3(/i,,//, - BjiUj) 

This confirms the former intuitive derivation of the perturbation SH to the escape probability 



as given in ( 40 ) . 



C Errata to the second order Compton collision term 

Although the second order expansion of the Compton collision term is present in many ref- 
erences, we find that often there are typos (actually we did not find one reference which is 
fully correct). For this reason, we prefer to summarize the corrections we found to what is 
reported in |32] and in [Tj)]. Let us start with Ref. [4"T] . 



42 



• In the third line of (2.9) in [37], there is an extra prime in the second factor of the last 
term, i.e. it should read ... — f^(p')[l + f^°\p)] ■■■■ 

• In the second line of (2.13) in |47j . the last term in the curly brackets should be multiplied 
by T e , so it should read {... - T e [f°\p') - f (0) (p)]d5(p - p')/dp'} 

• In (2.15) in [37], the first factor should not be present, i.e. the equation should read 
c { ? v (p,i?) = [f {1) (i?)-f m (P)} - 

Using the Compton collision term from [37], with the above corrections, we checked the 
final expression for C(p) given in [TO] by their eq. (4.42). There is one correction that needs 
to be made there: 

• If one looks at eq. (4.12) of [19j, some of the above terms have been corrected, but 
still in our opinion there are some typos remaining. Still, in eq. (4.42) only one typo is 
remaining. This is the last term, the Kompaneets one, which should be multiplied by 
m e /p 2 , giving the standard ^7^2^ •••• 

For computing the temperature fluctuations at second order, we find it useful to write the 
Collision term in a generic frame: 



(2) I — 
C 'i ^ - V 2) M-- V, - F (2) Yo (n)--F^(v) + 
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Here, following the notation of BMR1, we have expanded the distribution function as 

F = + + l -F^ , (140) 

where is the first order perturbation and is the second order one. The expansion in 
sperical harmonics of a quantity f(n) is given by: 



, 121 + 1 



fim = i-ir\l-^- I dnf(h)Y^(n) . (141) 
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